/** 湿空气的一些常数 */
export const airConstants = {
  /** 开氏温度；单位：K */
  T0: 273.15,
  /** 大气压力；单位KPa */
  Pressure: 101.325,
  /** 干空气气体常数；单位：J/(Kg·K) */
  Ra: 287.042,
  /** 水蒸气气体常数；单位：J/(Kg·K) */
  Rv: 461.524,
  /** 分子量比值；molecular weight ratio */
  MWR: 0.621945
} as const

/**
 * ### 根据干球温度计算饱和水蒸气分压力；
 * @reference ASHRAE. 2013 ASHRAE Handbook Fundamentals SI Edition [M].
 * @see https://www.buildenvi.com/zhuanti/abc/skq/150424a3
 * @param t 干球温度，单位：℃
 * @returns 饱和水蒸气分压力，单位：kPa
 */
export function vaporSaturatedPressure(t: number) {
  let result
  const c1 = -5.6745359e3,
    c2 = 6.3925247,
    c3 = -9.677843e-3,
    c4 = 6.2215701e-7,
    c5 = 2.0747825e-9,
    c6 = -9.484024e-13,
    c7 = 4.1635019,
    c8 = -5.8002206e3,
    c9 = 1.3914993,
    c10 = -4.8640239e-2,
    c11 = 4.1764768e-5,
    c12 = -1.4452093e-8,
    c13 = 6.5459673
  const Ta = t + airConstants.T0
  if (t < 0.01) {
    result =
      0.001 *
      Math.exp(
        c1 / Ta +
          c2 +
          c3 * Ta +
          c4 * Math.pow(Ta, 2) +
          c5 * Math.pow(Ta, 3) +
          c6 * Math.pow(Ta, 4) +
          c7 * Math.log(Ta)
      )
  } else {
    result =
      0.001 *
      Math.exp(
        c8 / Ta +
          c9 +
          c10 * Ta +
          c11 * Math.pow(Ta, 2) +
          c12 * Math.pow(Ta, 3) +
          c13 * Math.log(Ta)
      )
  }
  return result
}

/**
 * ### 湿空气水蒸气分压力计算；单位：kPa
 * @param t 干球温度，单位：℃
 * @param tw 湿球温度，单位：℃
 * @param atm 大气压力，单位：kPa
 * @returns 水蒸气分压力计算；单位：kPa
 */
export function vaporPressure(
  t: number,
  tw: number,
  atm = airConstants.Pressure
) {
  if (tw > t) {
    throw new RangeError('tw must be less than t')
  }
  let d: number
  const pvs = vaporSaturatedPressure(tw)
  const ds = moisture(pvs, atm)
  if (tw <= 0) {
    d =
      ((2830 - 0.24 * tw) * ds - 1.006 * (t - tw)) /
      (2830 + 1.86 * t - 2.1 * tw)
  } else {
    d =
      ((2501 - 2.326 * tw) * ds - 1.006 * (t - tw)) /
      (2501 + 1.86 * t - 4.186 * tw)
  }
  return atm * (d / (airConstants.MWR + d))
}

/**
 * ### 含湿量是指每千克干空气中所含的水蒸气量，单位为kg/kg（a）
 * #### 其中18.015 268是水蒸气的分子量，28.966是干空气的分子量。
 * #### 通过将这两个数值相除，得到的结果是0.621 945，这是一个经验性的系数，用于计算湿空气中的水蒸气分压力与干空气分压力的关系‌
 * @param pvs 饱和水蒸气分压力，单位：kPa
 * @param atm 大气压力，单位：kPa
 * @returns 含湿量，单位：Kg/kg
 */
export function moisture(pvs: number, atm: number) {
  if (pvs > atm) {
    throw new RangeError('pvs must be less than atm')
  }
  return (airConstants.MWR * pvs) / (atm - pvs)
}

/**
 * ### 根据 ! 露点温度 ! 计算含湿量；单位：g/Kg;
 * @param td 露点温度，单位：℃
 * @returns 含湿量，单位：g/kg
 */
export function moistureFromDewPoint(td: number) {
  const c1 = 9.164e-3,
    c2 = 1.446e-4,
    c3 = 1.741e-6,
    c4 = 5.195e-8
  return (
    3.703 +
    0.286 * td +
    c1 * td * td +
    c2 * Math.pow(td, 3) +
    c3 * Math.pow(td, 4) +
    c4 * Math.pow(td, 5)
  )
}

/**
 * ### 湿空气焓值计算；
 * #### 干空气的定压比热Cp·a=1.006 kJ/(kg·℃)，水蒸气的定压比热Cp·v=1.86 kJ/(kg·℃)；
 * #### 湿空气的比焓一般是以1kg干空气作为基数进行计算的，伴随着1kg干空气的还有d kg水蒸气，如果取0℃的干空气和0℃水的比焓为零；
 * #### 2501kJ/kg是每千克0℃的水变成0℃的水蒸气所需要的汽化潜热；
 * @param t 干球温度，单位：℃
 * @param d 含湿量，单位：Kg/kg
 * @returns 焓值，单位：kJ/kg
 */
export function enthalpy(t: number, d: number) {
  return 1.006 * t + d * (2501 + 1.86 * t)
}

/**
 * ### 湿空气密度计算；单位：Kg/m3；
 * #### 工程中密度取1.2kg/m3 已足够精确
 * @param t 干球温度，单位：℃
 * @param pv 水蒸气分压力，单位：Pa
 * @param atm 大气压力，单位：Pa
 * @returns 密度
 */
export function density(t: number, pv: number, atm: number) {
  if (pv > atm) {
    throw new RangeError('pv must be less than atm')
  }
  const ta = t + airConstants.T0
  return (0.003484 * atm - 0.00134 * pv) / ta
}

/**
 * ### 湿空气的露点温度计算；单位：℃
 * @param pv 水蒸气分压力，单位：Pa
 * @returns 露点温度，单位：℃
 */
export function dewPointTemperature(pv: number) {
  let res: number
  const y = Math.log(pv / 1000)
  if (pv > 622.082) {
    res =
      6.54 +
      14.526 * y +
      0.7389 * Math.pow(y, 2) +
      0.09486 * Math.pow(y, 3) +
      0.4569 * Math.pow(pv, 0.1984)
  } else {
    res = 6.09 + 12.608 * y + 0.4959 * Math.pow(y, 2)
  }
  return res
}

/**
 * ### 相对湿度计算；单位：%
 * @param t 干球温度，单位：℃
 * @param tw 湿球温度，单位：℃
 * @param atm 大气压力，单位：kPa
 * @returns 相对湿度，单位：%
 */
export function relativeHumidity(
  t: number,
  tw: number,
  atm = airConstants.Pressure
) {
  if (tw > t) {
    throw new RangeError('tw must be less than t')
  }
  if (t === tw) {
    return 100
  }
  return (100 * vaporPressure(t, tw, atm)) / vaporSaturatedPressure(t)
}

/**
 * ### 计算湿球温度
 * @param h 焓值，单位：kJ/kg
 * @param atm 大气压力，单位：kPa
 * @returns 湿球温度，单位：℃
 */
export function wetBulbTemperature(h: number, atm: number) {
  // 初始化低、高和湿球温度
  let low = -100,
    high = 100,
    tw = 0
  // 设置精度
  const acc = 5e-4
  // 循环1000次
  for (let i = 0; i < 1000; i++) {
    // 计算中间值
    const dt = (high - low) * 0.5
    tw = low + dt
    // 计算饱和蒸汽压力
    const pvs = vaporSaturatedPressure(tw)
    // 计算湿度
    const ds = moisture(pvs, atm)
    // 计算湿球温度的焓
    const hmid = enthalpy(tw, ds)
    // 判断焓的大小
    if (h - hmid > 0.0) {
      // 如果焓大于湿球温度的焓，则更新低值
      low = tw
    } else {
      // 否则更新高值
      high = tw
    }
    // 判断精度是否满足要求
    if (Math.abs(dt) < acc) {
      break
    }
  }
  // 返回湿球温度
  return tw
}
